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I. INTRODUCTION 

Clustering plays a relevant role in condensed matter physics, particularly in relation with 
transition phenomena, as it has been recognized since more than sixty years ago when Bijl, 
Band and Frenkel (among others)0"i introduced the concept of real clusters, in place of the 
mathematical ones that Mayer had considered in the virial expansion of imperfect gasesS. 
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The apphcation of statistical mechanics formalism to describe clustering in equilibrium 
classical systems is mainly due to Hill In Hill's theory the concept of cluster is directly 
related to that of connectivity: two particles belong to the same cluster if they are connected 
through a path of directly connected particles. Therefore a crucial point in the identification 
of the clusters is the definition of directly connected particles (a bonded pair). Whereas 
thermodynamic properties are not affected by the particular definition used to identify the 
clusters, clustering and percolation properties do are very sensitive to this choice insteadi. 

Hill's definition of a bonded pair is based on an energetic criterion: two particles are 
bonded, in a given configuration, if their relative kinetic energy is less than minus the pair 
potential energy!. Frequently, however, a geometrical criterion is used instead of the ener- 
getic one. Thus, Stillinger criterion states that two particles are bonded if they are separated 
by a distance shorter than certain value (i, which is called the connectivity distanceil. 

In the past few years, several papers have renewed interest in the energetic and dynamic 
aspects of bonding, particularly in relation to the properties of wateiS. Hydrogen-bond 
life times have been calculated and new connectivity criteria have been proposed for water 
molecules. These works suggest that a bonded pair must be defined as two water molecules 
that are in some appropriate geometrical arrangement at least during a time interval of the 
order of the estimated hydrogen-bond life times. Recently has been pointed out the necessity 
of including dynamical information in the definition of the clusters considered in nucleation 
studiesUl^. It is also expected that taking into consideration the dynamics of aggregates 
would improve the description of some phenomena as the conductor /insulator transitions in 
microemulsions0 . 

Hill's theory separates the Boltzmann factor e(ri, = exp[—(3v{ri, r2)] into bonded (f) 
and unbounded (*) terms: e(ri,r2) = e'''(ri,r2) + e*(ri,r2). In this paper we will consider 
only systems interacting via pairwise additive potentials and will denote with f(ri,r2) the 
pair potential. As usual jS = l/ksT with T the absolute temperature and ks the Boltz- 
mann constant. Since e^(ri,r2) represents the basic probability density that two particles 
at positions riand r2, respectively, are directly connected (bonded), this separation gets a 



diagrammatic expansion for the partition functions in terms of real clusters instead of the 
mathematical clusters of Mayer.. 

Fugacity and density expansions similar to those obtained by Mayer and MontrolS for 
the ordinary pair correlation function 5f(ri,r2) has been found, within Hill's formalism, by 
Coniglio and collaborators for the pair connectedness function g''{ri,r2). This function is 
proportional to the joint probability density of finding two particles connected (directly or 
indirectly) and at positions ri and r2, respectively^. Moreover, by collecting nodal and 
non-nodal diagrams in these expansions an Ornstein-Zernike like relationship is obtained. 
Consequently, integral equations for (7"l'(ri,r2) can be posed0. These equations, as well as 
computer simulations, have been applied together with Stillinger criterion to study clustering 
and percolation in several continuum systems. They include: the ideal gasEUll; simple hard 
sphereslli@; hard spheres with (sticky) adhesionlll'0; hard spheres with square wellsil'ii; 
hard spheres with Yukawa tailsB; charged hard sphereS; Lennard- Jones fluids^ and dipolar 
hard sphereJ^^^. 

Hill's work is the basis for all the present clustering theories in continuum systems. 
Nevertheless, it has some limitations. In one hand, as Hill pointed out in his paperi, the 
essence of his theory is the concept of bonded pair. This fact limits the kind of cluster we 
can study. In the other hand both, the energetic and the geometrical, criteria considered by 
Hill yield a poor description of a "real" bond. 

Let first consider the second problem. If we use an energetic criterion to define a bond, 
then the functions and e* must be momenta, as well as position, dependent. However, 
the functions and e* used by Hill are only position dependent, since the momenta have 
been averaged. A geometrical criterion is neither enough by itself to warrant the existence 
of a real bond. After a short time two neighbor particles can be much separated each other 
if the relative velocity is high. These features suggest that the bond definition should be 
modified. In Section II we introduce a bond criterion which incorporates a residence time 
into Stillinger criterion. We shall call "chemical clusters" the aggregates so defined. A 
connectivity theory, already presented in a previous paperS, that allows to handle these 



clusters will be reviewed in Subsection V-A. 

Let us now focus on the energetic criterion and take, for example, three particles which 
remain next one to the others. If the potential interaction between each pair is not strong 
enough to form a bond then, according to Hill's theory, we will have three clusters of one 
particle each. Nevertheless, since the three particles remain close each other, then the system 
can be considered as a single cluster of three particles "collectively or non-specifically pair- 
bonded" . We call this kind of aggregates "physical clusters" . 

With respect to the geometrical criterion, it could happen that it does not matter how 
the particles interact each other, if they are close enough, then they are bonded. Thus, 
according to the Stillinger's criterion, there will be three bonded pairs and the particles 
will form a single three-particle cluster. In this way the Stillinger's criterion partly avoid 
the limitation that the energetic criterion imposes in Hill's work. Due to this feature, 
Stillinger's criterion has given good results in nucleation studies!. However, it does not 
include dynamical information. As a consequence one could take as a single cluster two 
physical clusters which are close during a very short time. In Section HI we will give a 
definition for physical clusters which adds a dynamical restriction to the Stillinger criterion. 
This restriction is equivalent to require a minimum lifetime for the aggregates. In Subsection 
V-B we present a simplified theory to obtain connectivity functions and mean cluster sizes 
for physical clusters. 

It can be worth-mentioning that the difference between the chemical and physical clusters 
that we introduce is, in some sense, analogous to that existing between chemical (or specific) 
and physical (or non-specific) adsorption. 

We remark that in this work, the time is thought in the frame of the microscopical 
evolution of the system whereas its macroscopic (thermodynamic) state remains unchanged. 

Since the work of Coniglio and collaborators, the central object in clustering and perco- 
lation theory of continuum systems is the cluster pair correlation function (7^(ri,r2). It is 
the joint probability density of finding two particles that belong to the same cluster of kind 
7 (7 =Stillinger, chemical or physical) at positions riand r2, respectively. In particular, for 
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Stillinger clusters we have gstiu{ji,T^2) = r2)- From cluster correlation functions we 

can calculate the mean cluster size defined 



where N is the number of particles of the system and the integrations are over the whole 
system volume V. In Eq.(|lD we have explicitly indicated its dependence on the system den- 
sity p = N/V. Percolation of an ensemble of particles is concerned with the existence of 
clusters that become macroscopic in size. Then the percolation transition occurs at a critical 
percolation density pc which mathematically verifies 



The rest of this article is arranged in the following way. In Section II and III we introduce 
the chemical and physical clusters, respectively. Section IV will be devoted to discuss some 
results from molecular djTiamics simulations. In section V we present statistical mechanics 
theories to describe both kind of clusters. We conclude with some comments on potential 
applications of these new criteria. 



As it was mentioned, Stillinger's criterion does not describe a real bond in a completely 
satisfactory way since it does not include any dynamical information. For example, two 
particles which are separated by a distance shorter than d, at a given instant, are considered 
bonded even when their relative velocity is high enough to separate the particles in a short 
time. On the other hand, the analysis of the residence time of a molecule in the neighborhood 
of another one has produced important advances in the study of the hydrogen bond in 
wateifliffl. All those facts suggest to include a residence time to describe a real bond. In 
this direction we propose the following bond definition: 




(1) 



lim S'-y(p) = oo. 



(2) 



P^Pc 



II. CHEMICAL CLUSTERS 
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Two particles in the system are bonded at the instant t if they have been 
separated by a distance shorter than the connectivity distance d during the 



Then we define a cluster in the same way as Hill did it: 

A chemical cluster is a set of particles such that between any pair of them 



As we can see, the residence time modifies the bond definition but, once the definition 
of a bond is given, the same connectivity rule for the cluster identification is used. We call 
chemical clusters to this aggregates because the stable bonds are the essence of the cluster 
formation, like the covalent bond is the essence in the building of molecules^. In the limit 
case r = we recover Stillinger's definition. 

The numbers d and r are two parameters of the theory that must be obtained from 
physical considerations on the system and the phenomenon under study. For example, in 
nucleation studies for Lennard- Jones fluids, where the particles interact via the pair potential 



is generally adopted the value d = 1.5crE3. According with its role in the theory, one possible 
election of r could be done on the basis of considering the mean value of the relative velocity 
between two atoms in the fluid at temperature T: 



where m denotes the mass of the atoms. Thus we choose r as the necessary time for the 
particles travel, in average, the connectivity distance: 



To have an idea of magnitude orders, let consider argon at T = 300 K. We have a = 
3.41 A ,e = 1.65 x lO^^iJ and m = 6.62 x lO^^e j^g^ ^-^^^ ^ ^ a and r = 0.83 ps or. 



time interval [t — r, t], where r is the residence time. 



exists a path of bonded particles. 




(3) 





in reduced units, r* = ra ^Je/m = 0.39; d* = d/a = 1.5. 
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It should be remarked that this estimation for r is just an example. The residence 
time will depend on the system and the particular phenomenon to be studied and can be 
determined by very different causes. In section VI we briefly comment on the application of 
chemical clusters to the case of water in oil microemulsions. There we discuss in particular 
on the determination of the parameters. 

III. PHYSICAL CLUSTERS 

Here, we will consider how to identify physical clusters, i.e. those aggregates in which 
specific bonds between pairs of particles are not necessary for the clusters existence. As we 
have mentioned, this kind of clusters could be of interest in the simulation of nucleation 
process. The Stillinger's criterion and some other improved criteria are currently used in 
those studies!. 

We propose, in order to define physical clusters, a criterion in which the dynamics of the 
system is implicitly involved: 

A set of particles form a physical cluster in the instant t if all they remained 
connected during the whole time interval [t — r, t]. 

Here the term "connected" must be understood in the Hill-Stillinger sense, i.e. two 
particles are connected if between them there exists a path of directly connected particles. 
Two particles being directly connected if they are at a distance shorter than d. 

The previous definition is incomplete. It must be accompanied by the additional con- 
dition that those particles of the system which do not belong to a given cluster can not 
contribute to the connectivity of that cluster. Actually, the best way of stating the crite- 
rion is operationally. In this vein it is convenient to introduce the concept of fragmentation 
event. A fragmentation event occurs when a cluster breaks up into two or more connected 
subaggregates. 

Suppose we identify the system clusters at t — r according to the ( "instantaneous") 
Stillinger's criterion. Select one of these clusters. When the system evolves from t — t 

7 



to t, their particles eventually connect and disconnect several times. If we consider just 
those particles which had belonged to the selected Stillinger cluster at t — r and consider 
the successive fragmentation events then, after the period [t — T,t], the particles will lie in 
several subaggregates. Each subaggregate forms a physical cluster at time t. 

In Fig. 1 we show a simple example where the differences among Stillinger, chemical and 
physical clusters can be appreciated. At the initial time t — t we show the instantaneous 
StiUinger clusters (a). At the final time t we identify Stillinger (b), physical (c) and chemical 
clusters (d). We can see that physical clusters are smaller in size and longer in life than 
Stillinger ones, which are instantaneous. Thus, physical clusters would be more stable 
entities than Stillinger ones. However, physical clusters do not need pair-bonded particles 
to exist whereas these last are essential for defining chemical clusters. In Fig. 1, particles 
1 and 3 form a physical cluster of two particles, but they do not form a chemical cluster: 
although they are directly connected at t they were not at t — t. 

An estimation of the characteristic time r and the connectivity distance d is necessary 
for physical clusters as well as for chemical clusters. In our calculations on Lennard- Jones 
systems, we use the same value of d that for the chemical clusters. In order to select a value 
for T it is necessary to estimate the mean hfe time of an aggregate. Nevertheless, we are 
interested in comparing the results of the two different criteria. Thus, we will use the same 
value of T for both, chemical and physical criteria. 

IV. MOLECULAR DYNAMICS 

A. Chemical clusters 

In order to show the implications of using the residence time restriction in the bond def- 
inition we have simulated by molecular dynamics (MD) a system of particles which interact 
through a Lennard- Jones potential. Monte Carlo simulation can not be used because the 
trajectory of the particles is required. 
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We use a leap-frog algorithm with velocity correction to simulate a system at constant 
volume and temperaturell. To identify the chemical clusters we modify the standard MD 
code in such a way that the resulting algorithm is as follow: 

1) All initial conditions are set for the MD (t = to)- 

2) The NEIGHBORS routine is called. It generates a table with all direct connections 
-in the Stillinger sense- between the atoms. The positions of the particles at t = to are used 
in that operation. 

3) The MD algorithm is carried out with a very short time step At until a pair of atoms 
disconnect in the Stillinger sense {t = ti). 

4) The REDUCTION routine is called. It checks and deletes from the neighbors list the 
direct connections that have disappeared. The position of the particles at t = ti are used in 
this step. 

5) If ti — to < T go back to step 3, else go to step 6. 

6) All clusters are separated by using the remaining neighbors list. At this point the 
neighbors list only contain the direct connections which survived a time r from t = to, e.g. 
the list of chemical bonds. 

7) All interesting information is saved (coordinates, cluster list, etc.) 

8) Initial conditions are set from the last configuration. Return to step 2. 

Steps 1 to 8 are the loop for the chemical cluster identification corresponding to the final 
configuration at time t = ti. 

In step 2 NEIGHBORS routine makes a list with the same structure of the Verlet's 
neighbors listil, but a cut off equal to the connectivity distance d is used. 

In step 3 the leap-frog algorithm is carried out with a time step At short enough as to 
allow that no more than one disconnection occurs. We chose At to find a disconnection in 
each 3At intervals in average. 

In step 4 the neighbors list created in 2 is reduced by deleting the disappeared connec- 
tions. The new connections created in the system are discarded. Because in the last At 
interval more than one disconnection can occur, we check all connections in step 4 and any 
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disappeared connection is deleted from the neighbor list. 

After the time interval ti — to grows (by successive iterations) beyond r, a list of chemical 
bonds (as defined in section II) remains in the neighbors list. Any standard algorithm for 
cluster identification can be used to separate each cluster knowing the bonded pairJ^BS. 
We have used Stoddard's oneilii. 



B. Physical clusters 

For physical clusters we use a similar method to that applied in the case of chemical 
cluster. Nevertheless, the clustering steps are modified to achieve the desired definition. 
The complete algorithm is: 

1) All initial conditions are set for the MD {t = to). 

2) Stillinger clusters are tabulated by standard routinesH. The positions of the particles 
at t = to are used in that operation. 

3) The MD algorithm is carried out in very short time steps At until a pair of atoms 
disconnect in the Stillinger sense (t = ti). 

4) The DIVISION routine is called. It checks each cluster in the table made in step 2 
and breaks up those clusters which were fragmented as a result of the disconnection in step 
3. The cluster table is updated. 

5) If ti — to < r go back to step 3, else go to step 6. 

6) Save the interesting information. The cluster table contains all the set of particles 
which are in agreement with the physical cluster definition of section III. 

7) Initial conditions are set from the last configuration. Return to step 2. 

Steps 1 to 7 are the loop to identify the physical clusters for the final configuration at 
time t = ti. 

In step 4 a clustering count must be done separately over each initial cluster identified in 
2. In this way the connections formed between different clusters in the time interval [to,ti] 
are not included in the count. 



10 



C. Results 



In this section we show the results obtained by molecular dynamics simulation for a 
Lennard- Jones fluid with the pair potential given by Eq.(^. Our aim is not to exhaust the 
study of the various aspects of the clustering that we will consider but instead is simply 
to remark the qualitative differences they show for the three kinds of clusters: Stillinger, 
chemical and physical. 

We will express all quantities in reduced units: p* = Na^/V, T* = ksT/e, t* = 
ra^^ ^Js/m and d* = d/a for density, temperature, residence and life time, and connec- 
tivity distance respectively. We have used the leap frog algorithm with velocity corrections 
at half step to produce a molecular dynamics with constant number of particles A^, volume 
V and temperature tS. 500 particles were placed in a cubic box with periodic boundary 
conditions and 1000 configurations were generated after stabilization to calculate time aver- 
ages. We use a cut off equal to 2.1a for the pair potential interaction and a tail correction for 
the remaining interaction. Because we use the same number of particles in all simulations, 
scaling laws were not used to extrapolate in the limit oo. The time step is varied 

between At* = 0.00001 and At* = 0.0001 depending on temperature to get a connection or 
disconnection approximately every three time steps. The parameters for the cluster iden- 
tifications were chosen d* = 1.5 and r* = 0.5 for chemical and physical clusters. In order 
to compare the results with the standard Stillinger criterion we use the same connectivity 
distance d* = 1.5 for the three criteria. 

In Fig. 2 the correlation functions for Stillinger clusters gsuii{T 1,2) ■, chemical clusters 
9chemij'i,2) and physical clusters gphys{.fi,2) are shown. We have used free boundary condi- 
tions in the calculation of the correlation functions for clusters^. It means that neither real 
position nor nearest image of a particle but the image belonging to the proper aggregate 
through the replicas is used. 

If two particles are separated by a distance shorter than the connectivity distance, then 
they belong to the same Stillinger cluster with certainty. Thus, for ri^2 < d, the function 
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gstiuij'1^2) coincides with the ordinary radial distribution function g{r 1^2) as it can be appre- 
ciated in Fig. 2. The observed discontinuity at ri 2 = d for all cluster definitions is typical: 
the probability for two particles to be connected at ri^2 > d-, even for ri^2 d'^, depends 
on the presence of a intermediate third particle and thus the probability of belonging to 
the cluster notably disminishes. The function gstiuij'1,2) is long ranged because the density 
of the system considered in the figure coincides with the percolation density for Stillinger 
clusters at the given temperature. The mean Stillinger cluster size Sstm diverges at this 
density. 

As it is expected gchem{ri^2) and gphys{ri^2) are smaller than gstuiir 1^2) for all ri^2 because 
a dynamical restriction is required in addition to the geometrical one. Comparing gchemij'1,2) 
against gphys{ri^2) in Fig. 2, we see that gchem{ri,2) is above gphys{ri,2) for every ri_2. That 
is just what one expect when the same value for r is used in both, chemical and physical 
clusters definition, because a life time restriction on the aggregate is a condition weaker than 
imposing a residence time restriction on the bond. 

In Fig. 3 we present the system coexistence curve obtained by Panagiotopoulo£il together 
with the percolation curves for the various cluster definitions. These curves separate the 
phase diagram in two parts: percolated (high densities) and non-percolated (low densities); 
and intersect the coexistence curve. 

To determine if a given configuration is percolated we look for a cluster which span 
through the replicas. For a given temperature we increase the density until get 50% of 
percolated configurations. In that case we say that the percolation density was reached^. 
We observe that the percolation curve for physical clusters is slightly less steep than for the 
Stillinger ones. Nevertheless, it is at higher densities as is expected. 

In the other hand, the percolation curve for chemical clusters is at higher densities 
than for the physical ones and has a pronounced temperature dependence. We will discuss 
in section VI how this temperature dependence of chemical clusters makes them special 
candidates to describe the metal-insulator transition in water in oil microemulsions. 

Fig. 4 shows the size distribution function for Stillinger clusters nstiii{s), chemical clus- 
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ters nchem{s) and physical clusters nphys{s). These functions measure the number of clusters 
of size s per particle. For both, chemical and physical clusters, the number of small clusters 
is larger than for the Stillinger ones. In the displayed example, T* = 1.4 and p* = 0.155, 
the biggest chemical clusters are approximately of size 60, the biggest physical clusters are 
of size 200 and the biggest Stillinger clusters are of size 400. The life time restriction in the 
physical clusters avoid the unphysical big aggregates yielded by the Stillinger criterion. This 
fact suggest that physical clusters could be useful in nucleation studies of saturated vapors. 

In Fig. 5 the curves nstiu{s), nchem{s) and nphys{s) are shown at the corresponding 
percolation density for each cluster definition. At the percolation density we obtain the 



in a Lennard- Jones system z = 2.1 ± 0.1^. We found, averaging over all temperatures, 
zstiii = 2.05 ± 0.05, zchem = 2.14 ± 0.05, zp^ys = 2.30 ± 0.05. As we can see, for chemical 
clusters we have, within the statistical errors, a similar exponent that for the Stillinger ones, 
whereas for physical clusters we observe that the exponent value differs. We obtained the 
value of zphys for physical cluster with a longer life time, r* = 1.0, at T* = 1.4, for which 
percolation density is p* = 0.2 ± 0.01 and zphys = 2.41 ± 0.04. This show that the critical 
exponent z is very sensitive to the life time selected for physical clusters. 

Looking for structural information of the clusters we have calculated its moments of 
inertia. Surprisingly, we get elongated clusters, no matter which definition is used. Two 
of the principal moments of inertia (Ii and I2) are always similar between them and larger 
than the third one (I3). This result was explained by Laria and Vericatil when small clusters 
are under consideration. However, they got approximately spherical big clusters for dipolar 
hard spheres. In Fig. 6 the ratio (Ji + I2) / (2/3) is plotted as function of the cluster size for 
each cluster definition. Physical clusters are the more spherical ones, with a more near one 
ratio. Because big clusters are rare we average (Ji + 12)/{21^) for all clusters bigger than 20 



power 



lawffl: 



n(s) = As~\ 



(4) 
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particles size obtaining (Ji + l2)/{21^) = 2.52. This is far away from the spherical value. 

The origin of the difference between our result and the result obtained by Laria and 
Vericat lies in the fact that we use free boundary conditions in the calculation of the inertia 
matrix0. Unphysical structures are simulated when periodic boundary conditions are used 
as in the Laria and Vericat work. 

Recently, non-spherical clusters in nucleation studies have received special interest. An 
experimental result shows that spherical particles aggregate in non-spherical clusters at 
the liquid-solid transitionlll. However, spherical clusters are expected for the gas-liquid 
transition^. Our simulation results show that at liquid densities non-spherical clusters are 
the more frequent case. 

We give here a very simple explanation, based on steric considerations, for the existence 
of non-spherical clusters. For simplicity, suppose a quasi-circular two-dimensional cluster 
[see Fig. 7(a)]. If we add a particle in the accessible surface of the cluster a symmetry break 
down will occur. Now, the accessible surface is greater on the right side of the cluster [see 
Fig. 7(b)]. This increases the probability that a second particle added at random lands 
in the right side of the cluster. So, the cluster tends to grow in an elongated structure. 
Nevertheless, this process can not go on for ever. If the cluster is very elongated, the lateral 
accessible surface [top and bottom part in Fig. 7(c)] is very big. Then, the probability 
that a new particle added at random lands in the lateral region dominate. The mean shape 
of a cluster will balance the two competing effects. Suppose a cluster with the shape of a 
spherical caped cylindre (cylindre length L and caps radius R) (see Fig. 8). The balance 
of the mentioned competing effects means that the lateral surface of the cylinder has to be 
equal to the surface of the caps. So, we chose L = 2R. A solid revolution ellipsoid with 
semiaxis R and 2R has a moments of inertia ratio {I1 + I2) / (2/3) = 2.5, which is quite similar 
to the value obtained from our simulations for clusters bigger than 20 particles. 

Another structural parameter is the radius of gyration Rg. It is defined for a cluster of 
size s as@ 



14 



where r j is the position of the ith particle in the cluster and is the position of the centre 
of mass of the cluster 

1 

rem = - XI (6) 

^ i=l 

In Fig. 9 we present the radius of gyration as a function of the size of the cluster at the 
percolation density for each cluster definition. As usual Rg follows the power law^i 

Rg oc s^/^/ (7) 

where -D/ is the fractal dimension of the clusters. In our simulation that power law applies 
only for clusters bigger than 5 particles size. The fractal dimension is the same for all 
clusters definitions. The chemical and physical clusters are more compact than Stillinger 
clusters because they have a smaller radius of gyration. In the calculation of Rg we used 
free boundary conditions as for all the previous structural quantities. 

V. THEORY 
A. Chemical clusters 

In this section we will summarize the main results of the theory that we have developed 
elsewhere^ to describe the clustering and percolation for chemical clusters. 

For a system of classical particles interacting via a pair potential f (rj, r^) we define a 
density correlation function p(ri, r2, pi, P2) which is A^(A^ — 1) times the probability density 
of finding two particles at phase space points (ri, pi) and (r2, P2) respectively: 

p(ri,r2,Pi,P2)=^3j^^^^^^^^ (8) 



X 



N 2 N N 

n exp[-/3|^] n n ^M-^viv,. v,)]dv^'V-'. 
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Here h is the Planck constant and Q^iVyT) the canonical partition function of the system. 
Then, in the same spirit of Hill and Coniglio et. al.llll, we separate exp[— /3i;(rj, r^)] into 
connecting and blocking parts 

exp[-/3v{ri, rj)] = f{ri, rj, p^, p^) + f*{ri, rj, p^, p^) + 1 (9) 

Here /''"(rj, rj, pj, pj) represent the basic probability density that two particles at configu- 
ration (rj, rj. Pi, pj), are bonded. The shorthand notation /'''(rj, r^-, pj, pj) = f^j (7 = f, *) 
will be sometimes used. 

Substitution of Eq.(0) in Eq.(D yields 

A^(A^ - 1) r ^ / M 

p(ri, r2, pi, P2) = }^3NNiQ^^v,T) ^M-f^n^i, rs)] (10) 
where the sum is carried out over all possible arrangements of products of functions f-j and 

fk,l- 

It should be noted that the functions f-j and f*j can depend on momenta as well as on 
the positions of the two particles, but the sum of fjj and /j*j must be momenta independent 
in order that Eq.(|^) be satisfied. Except by this last condition, the functions fjj and f*j 
are otherwise arbitrary for thermodynamic purposes. Obviously we choose them in such a 
way that the desired definition of bonded particles is achieved: 



exp 





/t(r.,r„p.,p,) = r-^L-^7-^^^J'--^^' (11) 



otherwise 

f (r.'r,'p„p,)= (12) 
[exp[ fJv[r„rj)\ i otherwise 

where rjj(i) is the relative position of particles i and j at time t. We see that, in fact, 
Eq.(^ is satisfied by Eqs. (|Tl|) and (|1^). Explicitly, time is introduced here by taking the 



set {r^, p^} as initial conditions in t = and solving the equation of motion of the particles 
under their mutual interaction. 
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In order to exactly calculate rij{t) for any t we must solve a many body problem. An 
approximation can be obtained by reducing it to a two-body problem. This can be done by 
takingEi the effective force [— Vf ^-^-^(rj, rj)] composed of the direct force exerted by particle 
j on i [— Vf (fj, Tj)] plus the conditional average force exerted by the remaining particles in 
the system on the particle i. The last is an average over all configurations of the N — 2 
particles keeping i and j fixed. Thus 

\ k^ij I N-2 

It can be shownlll that this approximation corresponds to take f ^-^-^(rj, r^) = — ln[(7(rj, r^)]//?. 
This way rij{t) is obtained in terms of just the initial values rj, r^, pj and pj. 

Equation ([TT|) states that two particles i and j are bonded if they are separated by 
a distance shorter than d during a time interval longer than r. This coincides with the 
definition introduced in section II for chemical clusters. 

Each term in the integrand of Eq.(p!0|) can be represented as a diagram consisting of 
two white ei and e2-points, N — 2 black Cj-points and some f-j and /*^- connections except 

2 

between the white points. Here we take = exp [-/?|^]. White point are not integrated over 
whereas black points are integrated over their positions and momenta. All the machinery 
normally used to handle standard diagrams in classical liquid theoryS can be now extended 
to treat these new kind of diagrams. By following the Coniglio recipe to separate connecting 
and blocking parts in correlation functions i.e. g{ri, = g^{ri, r2, Pi, P2)+fl'*(i'i5 ^2, Pi, P2)j 
we obtain an Ornstein-Zernike like integral equation for g'^{ri, r2, pi, P2) 

g^ri, ra, pi, P2) = c^(ri, ra, pi, P2) (13) 

./p(.,P3)c,.,.,p.p,..(.,.,p.p,...p3. 

Here p(ri, pi)p(r2, P2)5'^(ri, ra, Pi, P2) is A^(A^ — 1) times the joint probability density 
of finding two particles at positions ri and ra with momenta pi and pa, respectively, and 
belonging to the same cluster, where the bonding criterion is given by Eqs.(|r^) and (P^; 
being 
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P(ri, Pi) = jy^Y / P(i"i' i"2, Pi, P2)c?r2rfp2. (14) 

The function c^(ri, r2, pi, P2) denotes the sum of all the non-nodal diagrams in the diagram- 
matic expansion of (7^(ri, r2, pi, P2)- We remember that a nodal diagram contains at least a 
black point through which all paths between the two white points pass. 
For homogeneous systems we have 

g^ir,, r2, Pi, P2) = c\r,, r^, Pi, P2) + (27rm/Cr)3/2 ^^^^ 

exp[-/3^]c"^(ri, rg, pi, P3)^^(r3, r2, Ps, P2)dr3dps. 

To get an integral equation from Eq.(^) is necessary a closure relation be- 
tween (^"''(ri, r2, pi, P2) and c"''(ri, r2, pi, P2)- We take the Percus-Yevick approximation 
(^(ri, r2) exp[/?t>(ri, r2)] = 1 + A^(ri, r2), where the function A^(ri,r2) is the sum of the 
nodal diagrams in the expansion of g{ri,r2). Separation into connecting and blocking 
parts i.e. 5f(ri, r2) = 5'"^(ri, r2, pi, P2) + g*{ri, r2, pi, P2) and A^(ri, r2) = A^"^(ri, r2, pi, P2) + 
A^*(ri,r2,Pi,P2), yields 

^^(ri, r2, pi, P2) = [/*(ri, r2, pi, P2) + l][^^(ri, r2, pi, P2) - c^(ri, r2, pi, P2)] (16) 
+ exp[/3v{ri, r2)]fi'(ri, r2)/^(ri, r2, pi, P2). 

Eq. (|13D closed by Eq. ([T6|) gives an integral equation for (7^(ri, r2, Pi, P2)- 

From the function (yf^(ri, r2, pi, P2) we define the pair correlation function for chemical 
clusters 

gchemiri, ra) = J p(ri, Pi)p(r2, P2)5'"^(ri, r2, Pi, P2)rfpirfp2. (17) 

It is the joint probability density of finding two particles that belong to the same chemical 
cluster at positions ri and r2, respectively. 

Then, the mean cluster size Schem is given by Eq.(|I|) and the percolation density is 
calculated from Eq.(^ with 7 = Chem: 



Schem — 1 + pj- J gChemijl^ r2)(iri(ir2, 



lim_ Schemip) = oo. 

B. Weak criterion for physical clusters 

In order to treat the physical clusters within a theoretical framework we will restrict 
ourself to a weak version of the definition introduced in section III. The weak criterion is 
estabhshed as follows 

Two particles of the system belong to the same physical cluster at the time t if 
at that time exists a path of directly connected particles which link them, and 
all particles forming the path, including the two particles under study, were 

connected in time t — t. 

Connected and directly connected are understood in the Stillinger sense. The difference 
with the strong criterion established in Section III is that in the weak version we require 
that the particles be connected just at the extremes of the time interval and not during the 
whole interval. The idea for the future is to use the results obtained from this definition 
within a path-integral like formalism in order to describe clusters according to the strong 
criterion. 

Consider the function if(ri,r2) introduced by Xu and Stell in their "percolation in 
probability" theory@. This function measures the basic conditional probability for two 
particles to be directly connected if they are at positions ri and r2, respectively. The role of 
this function within Hill formalism is that of an auxiliary function in the Boltzmann factor 
separation into bonded and unbounded terms: 

4,3 = H{Yi,Y.j)e{ri,rj) (18) 

Each choice of i?(rj, r^) corresponds to a different bond definition. For example, for the 
Stillinger's definition of directly connected particles we have 
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To define bonded particles according with the weak criterion of physical clusters, we use: 

H{r,,r,) = H{n^,;r) = { ^^^^^-^^ ' " , (19) 

[ ri,2 > d 

where /i^(ri 2; t) is the probability density of finding two particles connected at time t — t 
and in relative position ri 2 at time t. Clearly, this selection of H{r 12', t) considers as bonded 
any two particles which are separated by a distance shorter than d and have been connected 
(directly or not) at time t — r. We explicitly include r in the functions to stress the life time 
dependence. 

Once the direct bond is defined, a connectivity Ornstein-Zernike relation between the 
pair correlation function gphys{f'i,2', ^) and the direct correlation function for physical clusters 
cphys{ri,2', t) can be obtained in the standard form@@: 

9Physij'i,2] r) = Cphys{ri^2] '^)+ P J '^Physiri,3] '^)9Physij'3,2'^ r)dY^. (20) 

The function gphys{ji^2] t) is A^(A^ — 1) times the probability density of finding two particles 
in positions ri and r2 and belonging to the same physical cluster. 

In order to solve equation (|20D a closure relation is needed. For ri 2 < d, we see that 
9Physij'i.2] t) is exactly equal to ^h'^(vi^2', t), because the conditional probability density that 
two particles belong to the same physical cluster under the assumption they were connected 
at t — r is equal to one if they are separated by a distance shorter than d. For ri 2 > d we 
can use the Percus-Yevick connectivity closured. 

The calculation of /i^(ri 2;T) remains to be discussed. Consider the temporal correla- 
tion function introduced by Oppenheim and Bloom in the study of nuclear spin relaxation 
G(r'^ 2) 1^1,2; ^)^- This function is the probability density of finding two particles in relative 
position r'^ 2 = ^2~ ^'1 initial time and in relative position ri 2 = r2— riat time t later. 

It is given by 
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N N 

X 



EE^ [r'i,2 - (r.(0) - r,(0))] 5 [n.^ - (r,(t) - r,(t))] 



where is Dirac's delta function and (.) mean an average over all initial configurations 
in the canonical ensemble. In the case of an ideal gas, i.e. non-interacting particles of mass 
m, G{r[ 2, i'i,2 5 1) is exactly given by 



G'K2,ri,2,r) = -(^l exp 



4r2 



(22) 



For systems with non-null interaction some approximations for G{r[2,'''^i,2,t) have been 
reportecl&0. 
Now we define 

h\r,,,; t)= [ ^''f'\'^ G{r[^,, r,,„ T)dr[^,, (23) 

9vi,2) 

where gstiii{i^i2) is the pair correlation function for Stillinger clusters and g{r'i2) is the 
ordinary pair correlation function. The function h^{ri^2', t) is proportional to the probability 
density of finding two particles connected at time t — t and in relative position ri 2 at time 
t. Note that in the final time t the particles can be connected or not, but they must be 
in relative position ri 2, whereas in the initial time they must be connected whatever the 
relative position is. 

We have applied the previous formalism to an ideal gas. In this case (7(^1 2) = 1 and 
G'(ri,2,r'i,2,T^) is given by Eq.(|2|). 

The correlation function for Stillinger clusters gstiu{f^i.2) was obtained, for this system, 
by Chiew and Glandt in the Percus-Yevick approximation [c^(ri^2) = for ri^2 > 
We used this solution to obtain h^{ri^2',T) from equation (pB]). As we have mentioned, 
gphys{i^i,2] t) is known exactly for ri^2 < d- To close equation (|20| ) we take the approximation: 
Cphys{ri^2] t) = for ri 2 > d. By using the Baxter's factorization methodil and the Perram's 
algorithm!! we have solved the integral equation for gphys{ri,2'j t) iii ^1,2 > d. 
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In order to analyze the approximations we have introduced in our theoretical approach, 
we have also simulated by molecular dynamics the non-interacting system and have obtained 
the pair correlation function for physical clusters defined in the weak sense. The simulation 
algorithm to identify the weak physical clusters is slightly different from that considered in 
Subsection IV-B for the strong criterion. It is carried out as follow: 

1) Identify a Stillinger cluster at t — r. 

2) Move the particles of this cluster to the new positions in t without using periodic 
boundary conditions. In a ideal gas periodic boundary conditions are unphysical. 

3) Identify the Stillinger clusters in the new configuration for the particles selected in 1. 
These clusters are the physical clusters. 

4) Repeat from 1 until cover all initial Stillinger clusters. 

Fig. 10 shows the function gphys{i^i,2', t), obtained from our theory and simulation for the 
weak physical clusters at reduced density p* = pd? = 0.2 and various reduced life times r* = 
T / {d^/J3m) . When r* ^ we recover the Stillinger clusters gphysij'1,2] t = 0) = g sm 1^2)- 

In Fig. 11 we present the mean cluster size for the physical clusters as a function of r* 
for various values of p* obtained from Eq.([I|) with 7 = Phys: 



At low densities, p* < 0.3, the theory fits very well the simulation results. For higher 
densities the approximations used for gstiuiji,2) and cphys{ri^2] t) for ri_2 > d lead to poor 
results. 



In this work we have introduced two criteria for the cluster identification in continuum 
systems. Both criteria are based in the dynamics of the system. The first one uses a residence 
time in the definition of the bond between pairs of particles, the second one uses a life time 




VI. COMMENTS 



A. Summary and further developments 
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in the definition of an aggregate. Due to the quahtative features of the clusters yielded by 
the criteria we call them chemical and physical clusters, respectively. 

The main results that show the differences between chemical, physical and the traditional 
StiUinger clusters are collected below. 

Physical clusters: 

a) The percolation curve presents a slope similar to the corresponding curve for StiUinger 
clusters. However, this curve locates at higher densities than for the StiUinger case. 

b) The clusters are more spherical than the StiUinger ones for sizes smaller than 20 
particles. For clusters bigger than 20 particles size, physical and StiUinger clusters have the 
same shape and are always elongated. 

c) Big clusters are rare in comparison with StiUinger clusters. At percolation density 
a power law is followed by the cluster size distribution, but with a strongly r-dependent 
critical exponent. 

d) The fractal dimension of the clusters is the same that for StiUinger clusters but the 
physical ones are more compact. 

Chemical clusters: 

a) The percolation curve present a small slope as compared with the corresponding 
StiUinger curve. If the same value of r is used for physical and chemical clusters, then the 
percolation curve appears at higher densities in the chemical case. 

b) The sphericity of the chemical clusters is in between the physical and the StiUinger 
ones. 

c) At the percolation density a power law is followed by the cluster size distribution with 
an almost r-independent critical exponent. 

d) The fractal dimension of the clusters is the same that for StiUinger clusters but the 
chemical ones are as compact as the physical clusters. 

The definitions used in this work restricts us to use molecular dynamics simulation, 
because the trajectory of the particles must be known. To determine which particle belong 
to a given clusters at a given time t simulations need to be run between t — t and t in very 
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short time steps. This takes a large computing time. Thus, approximate theories are very 
much desirable to be developed. 

The solution of Eq. (p!3| ) for chemical clusters is a very difficult task. Tensors of rank 
9 must be used to storage correlation functions and a double convolution has to be solved. 
Recently, F. LadoS has developed a method to calculate the pair correlations for a polariz- 
able Lennard- Jones system which is based in expanding the relevant functions in orthogonal 
polynomials. The integral equation he has to solve is similar to Eq. ([T3| ) but with pi and 
P2 representing the induced dipolar moments of the particles 1 and 2 respectively. We are 
currently working in the application of Lado's code to the chemical clustering problem for 
Lennard- Jones particles. 

In the theory for physical clusters the main difficulty is to have a good approximation for 
G(r'^ 25 i"i,2, ^) and gstiuiji,2)- In one hand, the present theories to calculate the correlation 
function for Stillinger clusters do not provide it for densities higher than the percolation 
density. Because percolation densities are higher for physical than for Stillinger clusters 
the theory developed here is unable to give results for densities close to the "physical" 
percolation density. In other hand, approximations for (^(r'^ 2) 5^1,2, ^) in interacting systems 
are available just for the limiting cases t ~ and t ^ oo&^. 

Because the theory developed for physical clusters correspond to a weak criterion, an 
improved theory to describe the physical clusters introduced in section III is desirable. We 
hope that using a path integral approach in the interval \t — r, t] , where the weak criterion 
is used in each step of the path, the "strong" criterion could be reached. 

B. Potential applications 

Water in oil microemulsions are formed by mixing water, oil and a surfactant. Surfactant 
configuration favors the formation of inverted micelles with water inside and oil external to 
spherical structures^. When the micellar concentration is low, the system is a poor electrical 
conductor. Nevertheless, when the concentration of micelles increases beyond a critical value 



24 



a conducting behavior is present. This conductor- insulator transition has been explained 
by a continuum percolation mo del§i. When an aggregate of micelles spans the system, 
charge carriers hop from micelle to micelle in the aggregate and the system becomes a 
conductor. The process of passing a charged carrier from a given micelle to another one to 
which the former is directly connected requires some time. One can think of at least two 
ways of modeling this feature. In one of them the real pair attraction between micelles is 
exaggerated using, for example, sticky potentials to mimic the real intermicellar bonds. In 
the other one, we consider a pair potential closer to the real one and require the two micelles 
be connected during a time longer than r. This last point of view agrees with the idea of a 
chemical cluster. 

Most of the work on conductor/transition in microemulsions has followed the first ap- 
proach. Thus, the percolation curve for Stillinger clusters in a sticky hard spheres system 
was shown to fit the experimental conductor-insulator curve for water/decane/sodium-bis(2- 
ethylhexil)sulfosuccinateii. The curve has a peculiar behavior: it presents a small slope as 
compared with any other systemlilill. This behavior coincides with that presented by the 
chemical clusters in a Lennard- Jones system as we expect from the considerations above: 
the sticky spheres can be defined as bonded by just requiring a geometrical restriction in a 
given instant. As the interaction potential is infinitely deep at contact no residence time is 
necessary to impose in order to ensure a real bond. However, a continuum pair potential is a 
more realistic interaction for micelles and, in that case, the Stillinger criterion is not able to 
describe the conductor-insulator transition for water in oil microemulsions. In conclusion, 
the hopping time, e.g. the minimum time two micelles have to be close in order for a charge 
carrier to hop between them, seems to be necessary if realistic intermicellar pair potentials 
are going to be considered. The chemical clusters provide an accurate model to describe 
the conductor-insulator transition in this case. The values of d and r, for a chemical bond 
in Section II, can be used for water-in-oil microemulsions as the maximum distance and 
minimum time required to allow charge carriers to hop between two given micelles. 

As we mentioned, the concept of cluster is of fundamental importance in nucleation 
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studies on saturated vapors. These theories concern the growth of hquid-density-aggregates 
in a saturated vapor, which are the seed of a hquid phase. A comparative study of different 
clustering criteria was done recently by Senger et. al.li. The most frequently used criterion 
for clusters identification in nucleation studies is the Stillinger criterion and some improved 
versions of it. Nevertheless, this criterion can wrongly consider as part of a cluster some 
gas particles which are close to the cluster for a short time. In this way the cluster density 
will be lower than the density of the cluster obtained by considering that those particles do 
not belong to it. Also, the Stillinger criterion consider two physical clusters, which are close 
together for a very short time, as a single clusteiil. And a single cluster formed of two real 
clusters joined by only a point during just a very short period of time, will have a density 
of particles that is lower than the corresponding density inside each of the two separated 
clusters. 

Several modifications have been used to avoid this unphysical cluster countingiS. How- 
ever, the inclusion of dynamical information in the identification of clusters has been pro- 
posed as essential in selecting long-lived clusters in nucleation studieJ^l. Although we have 
studied systems which are in equilibrium, it is possible to use the previous clusters defini- 
tions in computer simulations of non-equilibrium systems also. We hope the physical clusters 
introduced in this work could contribute to a more appropriate identification of the relevant 
clusters in nucleation because their existence depends on a life time restriction. 
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Figures Captions 

Figure 1: Physical cluster identification, a) Position of the particles at t — r. The 
two Stillinger clusters are painted in different ways, b) Position of the particles at t. The 
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two new Stillinger clusters are painted in different ways, c) As in (b) but the four physical 
clusters are painted in different ways, d) As in (b) but the five chemical clusters are shown. 

Figure 2: Correlation function for the clusters at T* = 1.4 and p* = 0.155. Circles, 
triangles and diamonds correspond to Stillinger, chemical and physical clusters respectively. 
The doted line is the (ordinary) radial distribution function. The selected density is the 
percolation density for Stillinger clusters. 

Figure 3: Percolation curve for the Lennard- Jones fluid. Circles, triangles and dia- 
monds correspond to Stillinger, chemical and physical clusters respectively. Solid squares 
correspond to the coexistence curve^^ 

Figure 4: Log-Log plot of the cluster size distribution. Sohd, doted and dashed hne 
correspond to Stillinger, chemical and physical clusters respectively. T* — 1.4 and p* — 0.155 
in all cases. 

Figure 5: Log-Log plot of the cluster size distribution. Circles, triangles and diamonds 
correspond to Stillinger, chemical and physical clusters respectively. In each case the perco- 
lation density for T* = 1.4 was chosen. The z values were obtained by linear regression. 

Figure 6: Ratio between the largest (Ii and I2) and smaller (I3) moments of inertia 
as function of the cluster size. Circles, triangles and diamonds correspond to Stillinger, 
chemical and physical clusters respectively. The horizontal solid line is the averaged value 
over all clusters larger than 20 particles size. 

Figure 7: Two dimensional aggregates.. 

Figure 8: Average shape of a big cluster. 

Figure 9: Log-Log plot of the radius of gyration as a function of the cluster size at 
T* — 1.4 and percolation density. Circles, triangles and diamonds correspond to StiUinger, 
chemical and physical clusters respectively. The fractal dimension was obtained by linear 
regression. 

Figure 10: gphys{ri,2'i'^) vs. ri,2 for p* — 0.2, theory (solid line) and simulation (sym- 
bols). Diamonds, circles, triangles and squares correspond the r* = 0,0.1,0.316 and 1.0 
respectively. Note the scale change at ri,2/o? = 1. 
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Figure 11: Inverse of the mean cluster size for physical clusters as a function of r*, 
theory (solid line) and simulation (symbols). Squares, circles and triangles correspond to 
p* = 0.2, 0.3 and 0.5 respectively. 



31 





(a) 






1.'' •'tSt^^P 


(c) 


(d) 




s (cluster size) 




s (cluster size) 



5- 



O 



p*=0.155 
T*=1 .4 



CO 

C\J 4_ 



3- 



O 
A 
O 



o 

A 

O 



o 

A 

o 



o 



o 



o 



T" 
3 



T 
4 



"T 

5 



T" 

6 



T" 
7 



T 
8 



"T 
9 



10 



cluster size 



(a) (b) (c) 




1 




1 10 100 

s (cluster size) 



